{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.linalg\n",
    "import sympy as sy\n",
    "sy.init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> LU Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LU factorization rises because of computational efficiency, it mainly facilitates solving the system of linear equations, though the flops (number of floating operations) of LU is still higher than Guassian-Jordon. Nonetheless LU factorization is particularly handy if you are computing multiple solutions of $A x =b$. \n",
    "\n",
    "For intance, if you have a set of $\\{b_i,\\ i \\in \\mathbb Z\\}$ to substitute in the system, we only have to decompose $A$ once, but Guassian-Jordon requires the same amount of flops each time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> LU Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We aim to decompose a matrix $A$ into $L$ and $U$, which are lower and upper triangular matrix respectively. And $L$ must have $1$'s on its principal diagonal.\n",
    "\n",
    "$$\n",
    "A = LU\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason why $A$ can be decomposed is that $A$ can be converted into a upper triangular matrix after a series of row operations, i.e. \n",
    "\n",
    "$$E_p...E_2E_1  A = U$$\n",
    "\n",
    "Then \n",
    "\n",
    "$$A = (E_p...E_2E_1)^{-1}U = LU$$\n",
    "\n",
    "where \n",
    "\n",
    "$$L = (E_p...E_2E_1)^{-1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will hand calculate an example here:\n",
    "\n",
    "$$\n",
    "A =\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\3 & 4 & 6\\\\0 & 8 & 8\\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\-\\frac{1}{3} & 1 & 0\\\\0 & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\3 & 4 & 6\\\\0 & 8 & 8\\\\ \n",
    "\\end{bmatrix}}_{A}\n",
    "=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 8 & 8\\\\ \n",
    "\\end{bmatrix}}_{E_1A}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\0 &1  & 0\\\\0 & -\\frac{8}{3} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 8 & 8\\\\ \n",
    "\\end{bmatrix}}_{E_1A}\n",
    "=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 0 & -\\frac{8}{3}\\\\ \n",
    "\\end{bmatrix}}_{E_2E_1A=U}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "L^{-1} = E_2E_1 = \n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\0 &1  & 0\\\\0 & -\\frac{8}{3} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\-\\frac{1}{3} & 1 & 0\\\\0 & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1 & 0 & 0\\\\-\\frac{1}{3} & 1 & 0\\\\\\frac{8}{9} & -\\frac{8}{3} & 1\\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute $E_1^{-1}$ and $E_2^{-1}$, because\n",
    "$$\n",
    "L = E_1^{-1}E_2^{-1}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct augmented matrix for $E_2$ and $E_1$\n",
    "$$\n",
    "[E_1| I]=\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "-\\frac{1}{3} & 1 & 0 & 0  & 1 & 0\\\\\n",
    "0 & 0 & 1 & 0 & 0 & 1\n",
    "\\end{array}\\right]\n",
    "\\sim\n",
    "\\left[\\begin{array}{ccc|ccc}\n",
    "1 &0 &0 &1 & 0 & 0\\\\\n",
    "0& 1& 0& \\frac{1}{3} & 1 & 0\\\\\n",
    "0 & 0 & 1 &0 & 0 & 1\n",
    "\\end{array}\\right]\n",
    "=[I|E_1^{-1}]\n",
    "$$\n",
    "\n",
    "$$\n",
    "[E_2| I]=\n",
    "\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0  & 1 & 0\\\\\n",
    "0 & -\\frac{8}{3} & 1 & 0 & 0 & 1\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\sim\n",
    "\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0  & 1 & 0\\\\\n",
    "0 & 0 & 1 & 0 &\\frac{8}{3}& 1\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=[I|E_2^{-1}]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we get the inverse of $E_1$ and $E_2$:\n",
    "$$\n",
    "E_1^{-1} = \\left[\\begin{array}{ccc|ccc} 1& 0& 0 \\\\ \\frac{1}{3}& 1& 0\\\\ 0& 0 &1  \\end{array}\\right]\\\\\n",
    "E_2^{-1} = \\left[\\begin{array}{ccc|ccc} 1& 0& 0 \\\\ 1& 1& 0\\\\ 0& \\frac{8}{3} &1  \\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "L = E_1^{-1}E_2^{-1}= \n",
    "\\left[\\begin{array}{ccc|ccc} 1& 0& 0 \\\\ \\frac{1}{3}& 1& 0\\\\ 0& 0 &1  \\end{array}\\right]\n",
    "\\left[\\begin{array}{ccc|ccc} 1& 0& 0 \\\\ 1& 1& 0\\\\ 0& \\frac{8}{3} &1 \\end{array}\\right]\n",
    "=\n",
    "\\left[\\begin{array}{ccc|ccc} 1& 0& 0 \\\\ \\frac{4}{3}&1& 0\\\\ 0& \\frac{8}{3} &1 \\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We finished the LU decomposition\n",
    "$$\n",
    "A = LU =\n",
    "\\underbrace{\n",
    "\\left[\\begin{array}{ccc|ccc} 1& 0& 0 \\\\ \\frac{4}{3}&1& 0\\\\ 0& \\frac{8}{3} &1 \\end{array}\\right]}_{L}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\0 & 3 & 4\\\\0 & 0 & -\\frac{8}{3}\\\\ \n",
    "\\end{bmatrix}}_{E_2E_1A=U}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a look at SciPy results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9, 3, 6],\n",
       "       [3, 4, 6],\n",
       "       [0, 8, 8]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[9, 3, 6], [3, 4, 6], [0, 8, 8]]); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 0., 1.],\n",
       "       [0., 1., 0.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        , 0.        ],\n",
       "       [0.        , 1.        , 0.        ],\n",
       "       [0.33333333, 0.375     , 1.        ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[9., 3., 6.],\n",
       "       [0., 8., 8.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P, L, U = sp.linalg.lu(A)\n",
    "P\n",
    "L\n",
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9., 3., 6.],\n",
       "       [3., 4., 6.],\n",
       "       [0., 8., 8.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P@L@U # this is A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to notice the SciPy ```lu``` function does not only provide $L$ and $U$, it also has $P$ which is a permutation matrix. Theoretically, we don't need row switches to convert $A$ into $U$, but in some situations we must make row switches before, otherwise decomposition will fail to materialize.\n",
    "\n",
    "Thus Scipy uses $PLU$ decomposition to make the procedure always possible\n",
    "$$\n",
    "A  = PLU\n",
    "$$\n",
    "\n",
    "Actually $P = P^{-1}$, why? It is clearer in augmented matrix form, the inverse of row-switched elementary matrices are themselves.\n",
    "\n",
    "$$\n",
    "[P| I]=\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 0 & 1 & 0 & 1 & 0\\\\\n",
    "0 & 1 & 0 & 0 & 0 & 1\n",
    "\\end{array}\\right]\n",
    "\\sim\n",
    "\\left[\\begin{array}{ccc|ccc}  \n",
    "1 & 0 & 0 & 1 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0 & 0 & 1\\\\\n",
    "0 & 0 & 1 & 0 & 1 & 0\n",
    "\\end{array}\\right]=[I| P^{-1}]  =[P^{-1}|I]\n",
    "$$\n",
    "\n",
    "With these knowledge, what we are decomposing is actucally\n",
    "\n",
    "$$\n",
    "PA = LU\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\0 & -\\frac{3}{8} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\-\\frac{1}{3} & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 0 & 1\\\\0 & 1 & 0\\\\\n",
    "\\end{bmatrix}}_{E_0}\n",
    "\\begin{bmatrix}\n",
    "9 & 3 & 6\\\\3 & 4 & 6\\\\0 & 8 & 8\\\\\n",
    "\\end{bmatrix}\\\\\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\0 & -\\frac{3}{8} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\-\\frac{1}{3} & 0 & 1\\\\\n",
    "\\end{bmatrix}}_{E_1}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9& 3& 6\\\\0&8 &8\\\\ 3& 4& 6\n",
    "\\end{bmatrix}}_{E_0A}\\\\\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "1 & 0& 0\\\\0 & 1 & 0\\\\0 & -\\frac{3}{8} & 1\\\\\n",
    "\\end{bmatrix}}_{E_2}\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 &3 &6\\\\0 &8& 8 \\\\0& 3& 4\n",
    "\\end{bmatrix}}_{E_1E_0A}\\\\\n",
    "E_2E_1E_0A &=\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "9 &3 &6\\\\0 &8& 8 \\\\0& 0& -1\n",
    "\\end{bmatrix}}_{E_2E_1E_0A}=U\n",
    "\\end{align}\n",
    "$$\n",
    "Rearrange that we can see $PL$\n",
    "$$\n",
    "A = \\underbrace{E_0^{-1}}_{P} \\underbrace{(E_1^{-1}E_2^{-1})}_{L}U\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Solving Linear System by Using LU Factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the linear system：\n",
    "$$\n",
    "\\begin{align}\n",
    "3x_1-7x_2 -2x_3+2x_4&=-9\\\\\n",
    "-3x_1+5x_2 +x_3 &=5\\\\\n",
    "6x_1-4x_2 -5x_4&=7\\\\\n",
    "-9x_1+5x_2 -5x_3+12x_4&=11\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "In matrix form:\n",
    "$$\n",
    "\\underbrace{\\left[\\begin{array}{rrrr}\n",
    "3 & -7 & -2 & 2 \\\\\n",
    "-3 & 5 & 1 & 0 \\\\\n",
    "6 & -4 & 0 & -5 \\\\\n",
    "-9 & 5 & -5 & 12\n",
    "\\end{array}\\right]}_{A}\n",
    "\\left[\\begin{array}{r}\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "x_3 \\\\\n",
    "x_4\n",
    "\\end{array}\\right]\n",
    "=\n",
    "\\underbrace{\\left[\\begin{array}{r}\n",
    "-9 \\\\\n",
    "5 \\\\\n",
    "7 \\\\\n",
    "11\n",
    "\\end{array}\\right]}_{b}\n",
    "$$\n",
    "Perform $LU$ decomposition:\n",
    "\n",
    "$$\\underbrace{\\left[\\begin{array}{rrrr}\n",
    "3 & -7 & -2 & 2 \\\\\n",
    "-3 & 5 & 1 & 0 \\\\\n",
    "6 & -4 & 0 & -5 \\\\\n",
    "-9 & 5 & -5 & 12\n",
    "\\end{array}\\right]}_{A}\n",
    "=\\underbrace{\n",
    "\\left[\\begin{array}{rrrr}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "-1 & 1 & 0 & 0 \\\\\n",
    "2 & -5 & 1 & 0 \\\\\n",
    "-3 & 8 & 3 & 1\n",
    "\\end{array}\\right]}_{L}\\underbrace{\\left[\\begin{array}{rrrr}\n",
    "3 & -7 & -2 & 2 \\\\\n",
    "0 & -2 & -1 & 2 \\\\\n",
    "0 & 0 & -1 & 1 \\\\\n",
    "0 & 0 & 0 & -1\n",
    "\\end{array}\\right]}_{U}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace $A$ in the system, we have $L(Ux) = b$, solve this pair of equations\n",
    "\n",
    "$$\n",
    "Ly = b\\\\\n",
    "Ux = y\n",
    "$$\n",
    "\n",
    "Construct augmented matrix $[L|b]$ and $[U|y]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\underbrace{\n",
    "\\left[\\begin{array}{rrrr}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "-1 & 1 & 0 & 0 \\\\\n",
    "2 & -5 & 1 & 0 \\\\\n",
    "-3 & 8 & 3 & 1\n",
    "\\end{array}\\right]}_{L}\n",
    "\\underbrace{\\left[\\begin{array}{r}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "y_3 \\\\\n",
    "y_4\n",
    "\\end{array}\\right]}_{y}\n",
    "=\n",
    "\\underbrace{\\left[\\begin{array}{r}\n",
    "-9 \\\\\n",
    "5 \\\\\n",
    "7 \\\\\n",
    "11\n",
    "\\end{array}\\right]}_{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{array}{rrrr|r}\n",
    "1 & 0 & 0 & 0 & -9 \\\\\n",
    "-1 & 1 & 0 & 0 & 5\\\\\n",
    "2 & -5 & 1 & 0 & 7\\\\\n",
    "-3 & 8 & 3 & 1 & 11\n",
    "\\end{array}\\right]\\sim\n",
    "\\left[\\begin{array}{rrrr|r}\n",
    "1 & 0 & 0 & 0 & -9 \\\\\n",
    "0 & 1 & 0 & 0 & -4\\\\\n",
    "0 & 0 & 1 & 0 & 5\\\\\n",
    "0 & 0 & 0 & 1 & 1\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{array}{r}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "y_3 \\\\\n",
    "y_4\n",
    "\\end{array}\\right]=\n",
    "\\left[\\begin{array}{r}\n",
    "-9 \\\\\n",
    "-4 \\\\\n",
    "5 \\\\\n",
    "1\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we solve $Ux = y$, to show this in NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3.,  4., -6., -1.])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.array([[3, -7, -2, 2], [0, -2, -1, 2], [0, 0, -1, 1], [0, 0, 0,-1]])\n",
    "y = np.array([-9, -4, 5, 1])\n",
    "np.linalg.solve(U, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the process is right, this is the solution of the system and we can verify results by invoking ```np.linalg.solve()```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3.,  4., -6., -1.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[3, -7, -2, 2], [-3, 5, 1, 0], [6, -4, 0, -5], [-9, 5, -5, 12]])\n",
    "b = np.array([-9, 5, 7, 11])\n",
    "np.linalg.solve(A, b)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
